Viscoelasticity and Stokes-Einstein relation in repulsive and attractive colloidal glasses 
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We report a numerical investigation of the visco-elastic beliavior in models for steric repulsive 
and short-range attractive colloidal suspensions, along different paths in the attraction-strength vs 
packing fraction plane. More specifically, we study the behavior of the viscosity (and its frequency 
dependence) on approaching the repulsive glass, the attractive glass and in the re-entrant region 
where viscosity shows a non monotonic behavior on increasing attraction strength. On approaching 
the glass lines, the increase of the viscosity is consistent with a power-law divergence with the 
same exponent and critical packing fraction previously obtained for the divergence of the density 
fluctuations. Based on mode-coupling calculations, we associate the increase of the viscosity with 
specific contributions from different length scales. We also show that the results are independent on 
the microscopic dynamics by comparing newtonian and brownian simulations for the same model. 
Finally we evaluate the Stokes-Einstein relation approaching both glass transitions, finding a clear 
breakdown which is particularly strong for the case of the attractive glass. 



I. INTRODUCTION 

Understanding dynamic arrest in colloidal system is 
crucial in disparate technological applications (e.g. food 
industry yy, biomaterials|2j, painting). Development of 
basic science also requires a deeper understanding of the 
different routes and mechanisms leading to dynamic ar- 
rest (glasses and gels) [31 IH |SJ [B] . In this respect, model 
colloidal systems are playing a very important role due to 
their versatility. It is indeed possible to tailor the shape, 
size and structure of the colloidal particles making it pos- 
sible to design specific colloidal interaction potentials [7]. 
Furthermore, accurate experimental methods are now 
available for investigating the structure and the dynam- 
ics of colloids even at the single particle level [5]. Unex- 
pected novel behaviors regarding the glass transition have 
been theoretically predicted [U [THl EU [HJ [T^] and exper- 
imentally observed fT^, [151 El [III [HI [E] in the cases in 
which colloidal particles interact, beside the hard-core, 
via a short-range attractive interaction potentials (when 
the attraction range is about one tenth of the particle 
diameter or less). The predictions, based on applica- 
tion of the mode coupling theory for supercooled liq- 
uids (MCT) [20] suggest that the standard packing-driven 
hard-sphere glass transition transforms - discontinuously 
in some cases - into a novel-type of glass transition driven 
by the short-range attraction. The competition between 
the two different arrest mechanisms introduces slow- 
dynamics features which are not commonly observed in 
molecular and atomic systems. Experiments on solutions 
of (hard-sphere like) colloidal particles (either PMMA 
or polystyrene micronetwork spheres) in the presence of 
small non-adsorbing polymers [TH [151 HZ] have shown 
that there exists a window of polymer densities in which 
the mobility of the colloidal particles has a maximum for 



a finite value of polymer concentration. Moreover, for 
small and large polymer concentrations, the strength of 
the Qf-relaxation (the non-ergodicity parameter) is found 
to be very different, suggesting that indeed the visco- 
elastic response of the repulsive and attractive glass will 
also be significantly different. Molecular dynamics sim- 
ulations of short-ranged models [2 1, 22 , 23 , 2Ai have con- 
firmed the picture resulting from the theoretical predic- 
tions and validated by the experiments. A recent review 
can help summarizing the experimental and numerical 
studies in short-range attractive colloids [5J. 

The numerical results have been so far mostly limited 
to the study of self and collective properties of the density 
fluctuations. Despite the strong link with experiments 
and the relevance to industrial applications, the numeri- 
cal evaluation of the viscosity, 77, and viscoelastic proper- 
ties fi{ijj) have lagged behind, since significant computa- 
tional effort is requested for accurate calculation of i){uj), 
even more for states close to dynamical arrest. Experi- 
mentally, measurements of ry close to the repulsive hard- 
sphere glass transition show an apparent divergence, but 
there is no consensus on the functional form describing 
such increase [m [IB] . For colloidal gels, a power law di- 
vergence has been reported in connection to the gel tran- 
sition [27 . Theoretically, MCT predicts an asymptotic 
power law divergence, with identical exponent, of all dy- 
namical quantities with the distance from the transition, 
and hence r?, the time scale of the density fluctuations 
T and the inverse of the self diffusion coefficient \/Dq 
should diverge with the same critical parameters. 

In this article, we attempt a characterization of the 
viscoelastic properties of two different short-range at- 
tractive potentials (a polydisperse Asakura-Osawa and 
a square-well) along three different paths in the attrac- 
tion strength-packing fraction plane, which allow us to 
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access both the repulsion driven and attraction driven 
glass transitions with both systems. We show the di- 
vergence of the viscosity, as well as the diffusion coeffi- 
cient or structural relaxation time, as the repulsive and 
attractive glasses are approached. At high density, the 
isochoric path shows the reentrant glass; the viscosity 
increases about three orders of magnitude upon either 
increasing or decreasing the strength of attraction. 

The article is organized as follow: in Sec. II we intro- 
duce the numerical models and describe the methods to 
calculate the viscosity. In Sec. Ill we describe the paths 
investigated and provide some background information 
on the behavior of the diffusion and collective density 
fluctuations along these paths. In Sec. IV we discuss the 
observed behavior of the viscosity on approaching the re- 
pulsive and the attractive glass lines. In Sec. V, guided by 
theoretical MCT predictions for the viscosity, we provide 
evidence that the visco-elastic behavior close to the two 
different glass lines is controlled by density fluctuation 
of different wavelength. Finally in Sec. VI we report a 
study of the density and attraction strength dependence 
of the Stoke-Einstein relation. 



II. NUMERICAL SIMULATIONS 



A. Model A: Square well and Hard Sphere Binary 
Mixture 



We perform Molecular Dynamics (MD) simulations of 
a 50:50 binary mixture of 700 particles of mass m with 
diameters (taa = 1-2 and aBB = 1 (setting the unit of 
length). The particles interact through a hard core re- 
pulsion complemented by a narrow square well (SW) pair 
potential. The hard core repulsion for the AB interac- 
tion occurs at a distance aAB = (ctaa + <^bb)/'^- The 
SW potential is. 



Vsw{r) 




r < (Ji 
r > a,. 



(1) 



A,; 



where r is the distance between particles of types — 
A, B, the depth of the well uq is set to 1 and the widths 
Aij are such that Aij/{aij + Aij) = 0.03. Tempera- 
ture T is measured in units of uq {Ub — 1), the attrac- 
tion strength T = 1/T, time t in aBBim/uoY^^ ■ The 
use of a binary mixture allows us to suppress crystalliza- 
tion at high packing fraction = (paO'a + Pb<^b) ' ^Z^' 
where pi = Ni/L^, L being the box size and Ni the num- 
ber of particles for each species. The system undergoes 
phase separation into a gas and a liquid for large attrac- 
tion strength in a wide range of packing fractions |28| : 
the critical point is located roughly at Fc « 3.33 and 
(j)c ~ 0.27 (the latter is estimated from the Noro-Frenkel 
scaling [29] invariance close to the Baxter limit |30j). Pre- 
vious studies [231 HSl IS] of the same model allowed us to 
locate the dynamical arrest line and the spinodal curve. 



The 'numerical' glass line was determined by extrapola- 
tion via a power-law fitting of the normalized diffusion 
coefficient D/Dq, i.e. D/Dq ~ (0 — (j)g)'^ [31], where 
Dq = F^/^ . This study was complemented by the calcu- 
lation of the MCT glass lines for the same model. Hence, 
a bilinear transformation of and T was used to to su- 
perimpose the theoretical onto the numerical glass line. 

We also study, as discussed below, the same 
50:50 binary mixture of 700 particles, with the same 
<^AA, f^BB, o'AB above, but interacting simply as hard 
spheres, for which the potential reads. 



Vnsir) 







r < (Ji . 



(2) 



For Newtonian dynamics (ND) simulations, we used a 
standard event-driven (ED) algorithm [32]. We also per- 
form Brownian Dynamics (BD) simulations of the same 
model, to ensure the independence of the viscoelastic cal- 
culations on the microscopic dynamics. For BD simula- 
tions we exploit a recently developed BD algorithm, 
which we shortly describe below. For a more extensive 
discussion we invite the reader to consult Ref. [M]. 

If the position Langcvin equation is considered, i.e.: 



Do 



f,(t)+f,(t), 



(3) 



where Vi(t) is the position of particle i, Dq is the short- 
time (bare) diffusion coefficient, ii{t) is the total force 

acting on the particle, ri{t) a random thermal noise sat- 
isfying < ri{t) ■ ri{0) >= 6DoS{t). The BD integration 
scheme of Eq. [3] can be schematized as follow: 

(i) every t„ = nAt (n integer) extract velocities Vi 
according to a Maxwellian distribution of variance 
^/kBT/m■, 

(ii) evolve the system between t„ and tn + At according 
to the laws of ballistic motion (performing standard 
ED molecular dynamics). 

In other words, Gaussian particle displacements Afl — 
vlAt are extracted according to (Ar^^) — QDoAt and be- 
tween two velocities extractions, standard ED dynamics 
is applied. 

The present binary mixture model allows us to study 
the viscoelastic properties within the reentrant liquid re- 
gion, enclosed by the nearby attractive and repulsive 
glass transitions. On the other hand, due to phase sep- 
aration, it does not allow us to approach the attractive 
glass line at moderate density. Hence we will study Vhs 
for varying <j) (Path lA in Fig. [T]) and Vsw at fixed 
(/) = 0.58 on varying T (Path 3 in Fig.[l]). 

B. Model B: Asakura-Oosawa Polydisperse System 

We also study an interaction potential based on the 
Asakura-Oosawa model to make a direct link with ex- 
periments in colloid-polymer mixtures. A polydisperse 
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system, comprised of 1000 particles, is simulated with 
the standard velocity Verlet algorithm for Newtonian Dy- 
namics in the canonical ensemble, which requires a con- 
tinuous differentiable potential. To this end, a soft core 
was used instead of the hard core in Model A: 



singularity. We will use this system to approach the re- 
pulsive glass with increasing (j)c at (j)p = 0, hence using 
simply Vsc (Path IB in Fig. [l]), as well as to study the 
attractive glass at moderate density (j>c = 0.40 (Path 2 in 
Fig.[l| by using Vtot- 



Vscir) = {a,, It) 



36 



(4) 



C. Computation of viscosity 



where CTy — (oi ~\- aj)/2, with Ui the diameter of parti- 
cle i. Diameters where distributed according to the flat 
distribution [a — 5^a + S\ with a the mean diameter and 
5 = 0.1(7. The short-range attraction between particles 
is given by the Asakura-Oosawa model for polydisperse 
systems: 



(^ + lf-|(^ + l)V^ 



3^ . .2 



(5) 



for (Ti2 < r < <Ji2 + ^i^d for larger distances; 
Vi ~ (^il^i V — (^1 + ?72)/2, and (fip is the volume frac- 
tion of the polymer. The range of the interaction, ^, is 
the polymer size, and its strength is proportional to (jjp, 
the concentration of ideal polymers. To ensure that the 
interaction potential Vsc + Vao has its minimum at (T12, 
the Asakura-Oosawa potential is connected analytically 
to a parabola at ai2 -1-^/10 [35]. For average particles, 
(Ti = (72 = cr, the attraction strength of the Asakura- 
Oosawa potential is given by Vmin = —kBT(j)p{?i/2ri+ 1), 
which for ^ = 0.1, is Vmin — —l6kBT(j>p. 

Because the attractive glass transition occurs inside 
the liquid-gas spinodal, it cannot be accessed directly 
from the fluid with this potential. Thus, we have added 
a long range repulsive barrier to the interaction potential 
that destabilizes a macroscopic separation into two fluid 
phases. The barrier is given by: 



The shear viscosity t] is given by the Green-Kubo re- 
lation: 







(8) 

which expresses i] as the integral of the correlation func- 
tion of the non-diagonal terms of the microscopic stress 
tensor, cr"'' = YlfLi mViaVi^ - E^^j ^'(^'i)' 
where V is the volume of the simulation box, Via is the 
a-th component of the velocity of particle i, and V is the 
derivative of the total potential. (...) indicates an average 
over initial conditions. However, from the computational 
point of view it is more convenient to use the Einstein 
relation, 

^ = lim,y(i) = A lim j{AAm, (9) 

where AA{t) is the integral from to t of the three off- 
diagonal terms of the stress tensor. 



AA{t) = A{s + t)~- A{s) 



s+t 



Y,a'^P{s')ds' (10) 



Using Eq|9]is analogous to the calculation of the diffusion 
coefficient as the long time slope of the mean squared 
displacement. 

For discontinuous potentials (hard cores or square 
wells), equation [9] can still be used|31| despite the im- 
pulsive character of the interactions. In this case. 



Vhar{r) = ksT ■ 



r — ri 



(6) 



for Tq < r < r2 and zero otherwise, with ri = (r2 + ro) /2. 
The limits of the barrier were set to tq = cri2 -I- ^, and 
r2 = 2a, and its height is Ifc^T. The barrier raises the 
energy of a dense phase, so that liquid-gas separation is 
suppressed. The resulting total interaction. 



Vtotir) = Vscit) + VAoir) + Har(0 



(7) 



N 

[AA{t)]HS.SW = ^ [("^ ViaVi[))Tt + 

collisions 0.^(3 i — 1 

m{xka - a^;a)(Wfe^*'"' - v'il}'^'"''')] (11) 

where is the time elapsed from the previous collision, 
k and I are the two colliding particles, Xka is the position 



of particle k in direction a, and {v'^j-^'^^ 



V 



before 
IB 



) is the 



momentum change in direction /3 of particle k due to 
the collision with particle I. We have not attempted to 
numerically recover Caa{t) from AA{t). 



is analytical everywhere and allows straightforward inte- 
gration of the equations of motion. 

This model allows us to study the viscoelastic proper- 
ties of the fluid close to the attraction driven glass tran- 
sition at moderate density, i.e. far from the high order 



D. Units 

For both studied models we report states in the pack- 
ing fraction vs. attraction strength plane (0c — P)- For 
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4'° 


b 


It 


7d 


Model A: Vhs 


0.584 


0.51 


2.75 


2.17 


Model B: Vsc 


0.594 


0.53 


2.72 


2.02 



TABLE I: Glass transition point von Schweidler exponent 
b, and divergence exponents of the characteristic time of the 
decay of density fluctuations 7,- and of the diffusion coeffi- 
cient 7d for models A and B in the absence of attraction, i.e. 
respectively Vhs and Vsc, along path 1. 

Model A, the attraction strength is given by the inverse 
temperature (for HS temperature is irrelevant and is set 
equal to 1), whereas for Model B, F = —Vmin- Distances 
are measured using gb b for model A and the mean diam- 
eter, a for model B, while the particle mass, m, is always 
set to one. The stress correlation function is measured 
in units of fc^T/cr'^, and time in units of (a^m/kBTY^^ . 
The viscosity is thus given in (mkBT)^^^ /a^ . For the 
integration of the equations of motion in model B, the 
time step was set to St — 0.0025/-\/3. 

III. DESCRIPTION OF PATHS, TRANSITION, 
FITS, EXPONENTS 

Using the models presented above, we numerically 
study the following paths schematized in Fig. [T] 

Path 1: The zero-attraction case for both models, i.e. 
the hard- and the soft sphere models. The two models 
are not identical along this path because (i) the Asakura- 
Oosawa model has a soft repulsion (although the r^'^^- 
core is quite hard and no important effects are expected 
[37]) and more importantly ii) the size distributions are 
different: bimodal in model A vs. continuous in model 
B. Model B has been studied previously along this path 
monitoring the self-diffusion and the density correlation 
functions [38]. The glass transition points and the expo- 
nents controlling the power-law divergence of the struc- 
tural relaxation time scale, 7,-, and the diffusion coeffi- 
cient, 7d, as well as the von Schweidler exponent b (which 
provides a measure of the slow-decay of the density corre- 
lation function), are shown in Table |l] for both systems. 
The difference in the critical packing fractions can be 
attributed to the different size distributions of the two 
models. The exponents and on the other hand, 
are very similar in both models. 

Path 2: Approaching the attractive glass. This path is 
studied with model B, for which the liquid-gas transition 
is destabilized and the glass transition can be approached 
from the fluid. This path has been studied previously 
monitoring the density correlation functions |35[ 139] and 
the viscosity [40,^, and the glass transition is found for 
r*^ = 9.099; the associated von Schweidler and critical 
exponents are given in Table [IT] 

Path 3: The reentrant region and the approach to the 
attractive glass. This path is studied with model A, at 
(f>c = 0.58, a value well within the reentrant region j23]. 
The corresponding parameters for this path are provided 







b 7t 


7d 


Path 2: 


Vtot 


9.099 


0.37 3.23 


1.23 


Path 3: 


Vsw 


3.56 


0.33 3.75 


2.2 



TABLE II: Glass transition point P'', von Schweidler expo- 
nent b, and divergence exponents 7^ and 70 for models A 
and B in the presence of attraction, i.e. Vsw and Vtot, along 
respectively path 3 and 2. 

in Table [IT] At large temperature, the glass transition is 
approached but not reached because the studied packing 
fraction is close, but smaller than cf)'^ for Vhs, i-e. the 
path is parallel to the repulsive glass line in the limit 
T ^ 00. 

Note that, as predicted from MCT, the attractive glass 
shows lower von Schweidler exponents than the repulsive 
glass, for both paths and models, while 7^ is larger. This 
implies that the divergence of the time scale for struc- 
tural relaxation is more abrupt. For the square well mix- 
ture, quantitative results from simulations and MCT are 
available |31j. predicting the transition point at (f> — 0.58 
for Y'-^''^'-^'^ ~ 3.70, in quite good agreement with that 
estimated from the fits F'' ~ 3.56. For path 2 a quanti- 
tative comparison with MCT has been also recently per- 
formed |41j . showing that the driving mechanism for the 
slowing down observed in the simulation is driven by the 
short-range attractions (large-g modes of S{q)). 



paths 1A-B 




Colloid packing fraction, 

FIG. 1: Schematic phase diagram showing the attraction and 
repulsion driven glasses and the three paths followed in this 
work. Note that path 1 (infinite temperature limit) is studied 
within both models. The inset shows the three paths in the 
temperature-packing fraction representation. 



IV. VISCOSITY RESULTS 

In this section we study the viscosity along the three 
paths described above. 
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A. Hard and soft spheres: Paths lA and IB 




t 



FIG. 2: Upper panel: Stress correlation function Ccrcr{t) for 
Vac- The thin lines are empirical fittings to describe the data 
(see section |V] for details). Lower panel: Full lines are (3 < 
{AA{t))^ >/6Vt (from the Einstein relation Eq. ^ for aU 
studied </>c. For two specific values of cpc {4>c ~ 0.57 and (/)c = 
0.40) we also show ri{t) obtained using a direct integration of 
Ccrcr{t) (symbols), and integration of the fitting curves (dashed 
thick). Note that while ri{t) and P{AA{t)f/6Vt have the 
same long-time value, their time dependence is different. 

In Figure [2] we present, along path IB, the stress cor- 
relation function for Vsc at different concentrations (up- 
per panel) , and the integral of the squared non-diagonal 
terms of the stress tensor (lower panel). The correla- 
tion functions have been averaged over 5000 independent 
calculations. Note the progressive development of a two- 
step decay in C'a-ait) as the concentration increases and 
the glass transition is approached, with the second (struc- 
tural) decay of Caa{t) moving to longer and longer times. 
This implies that stress relaxes slower and slower, or 
equivalently that the system increases its ability to store 
the stress; i.e. the system becomes viscoelastic. Addi- 
tionally, it can be observed that Caa{Q) grows close to the 
transition. Both effects are responsible for the increase 
of the viscosity upon increasing the packing fraction, but 
the increase in the time scale is the one providing the 



leading contribution to the integral (see Eq. [8| . 

The integral of the stress correlation function is very 
noisy, and the numerical evaluation of the viscosity is 
more accurate if calculated using the Einstein relation 
(Eq. [9]), as shown in the lower panel of Fig. [2j For com- 
parison, the integral of the functional form used to de- 
scribe Caa{t) (see below) is also included for two state 
points. Note that all three quantities show the same long- 
time limit, i.e. the viscosity does not depend on the way 
it is calculated. At intermediate times, the integral of 
Ca-cr{t) and its fitting are in perfect agreement, but the 
integral of the fitted function is less noisy. Thus, we will 
calculate viscosities using the Einstein relation in Eq. |9] 

The viscosity, as given by the long-time plateau, grows 
with increasing particle density, as shown in Fig. [3] This 
increase is consistent with a power-law, diverging at the 
transition point estimated from the structural relaxation 
time and from the diffusion coefficient, (p^ = 0.594 [4Qj. 
The exponent for this power-law 7^ = 2.74 is similar to 
7t but different from 7/5, as reported in Table [l] 
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FIG. 3: Viscosity of soft (full black circles) and hard (empty 
red circles) spheres as a function of particle packing fraction, 
approaching the glass transition. Lines are power law fits 
to points with cj) > 0.50. The values of the critical packing 
fraction have been fixed to the previously determined values 
(see Table |l]), i.e. (j)^ = 0.594 and 0° = 0.584 for soft and 
hard sphere respectively. The corresponding fitting exponents 
7^ are 2.74 and 2.9. 

For hard spheres, path lA, we only show the integrated 
squared non-diagonal terms — obtained from Eq |ll| — in 
Fig. |4] These results are obtained averaging over 20 inde- 
pendent starting configurations and over time for a min- 
imum of 70rQ,, where is the density relaxation time at 
the wavelength corresponding to the nearest-neighbour 
peak. The behaviour of the curves is very similar to that 
shown above for model B, and the viscosity, also shown in 
Fig. |3] increases as the glass transition is approached. A 
power-law divergence with exponent 7^ ~ 2.9 is observed 
for the viscosity, with transition point at 0^ = 0.584, 
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slightly lower than for Vsc- The value of the exponent 
is, again, in good agreement with 7t- but quite different 
from 7£i. 




t 



FIG. 4: 13 < (Ay4(t))2 > /6Vt (with /3 = 1) for hard spheres, 
along path lA. 



^500, ■ ■ 




B. Attractive glass: Path 2 



In this section, we analyse the viscoelastic behaviour 
close to the attractive glass. As discussed above, for this 
purpose we use model B for which the liquid-gas separa- 
tion is suppressed by the presence of the added repulsive 
barrier, allowing for the study of low density (0c — 0.40) 
in a homogeneous system. In Fig. |5] we present again 
the stress correlation functions and the calculation of the 
viscosity by integrating the squared stress tensor non- 
diagonal terms. The attraction between particles induces 
a minimum after the short time (microscopic) relaxation, 
which introduces a negative correlation at intermediate 
attraction strengths. The origin of this minimum is sim- 
ilar to that in the velocity auto-correlation function, al- 
though here it is caused by stretching and rebound of 
the bonds. At high attraction strength, the correlation 
is positive again at all times, and after the minimum, 
Cacr {t) shows the development of a two-step decay and a 
large increase of the value at zero time Co-ct(O), similarly 
to the phenomenology observed for the repulsive glass. 
This indicates that the system is becoming solid-like. 

{{AA{t)y), shown in the lower panel of Fig. [s] grows 
dramatically upon increasing the attraction strength. 
The long time limit value, rj, is shown in Fig. |6] as a 
function of attraction strength. The data can be fitted 
using a power law divergence as a function of the dis- 
tance from the transition, F — F*^, where F*^ is reported 
in Table |TTj The exponent 7,, = 3.16 is again in good 
agreement with 7^-. 



FIG. 5: Stress correlation function Caaii) (upper panel) and 
P < (Ay4(t))^ > l&Vt (lower panel) for different state points 
along the isochore (j)c = 0.40. The thin lines in the upper 
panel represent empirical fittings to Ccra{t), eq. (see section 
|V]for details). 



C. Reentrance region: Path 3 

As discussed above, path 3 is a high density isochoric 
path, where the attractive and repulsive glass lines are 
about to merge. Varying the attraction strength, the sys- 
tem can be studied in states close to the repulsive or to 
the attractive glass. This path is studied only with sys- 
tem A, because the short interaction range of the studied 
SW opens up a large fluid region between the two glasses. 
Fig.|7]shows ((AA(t))2)/t calculated using Eql 
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_ _ The 

corresponding viscosity is reported in Fig. |6]as a function 
of F. As expected in this region, the viscosity increases 
both at low temperature, due to the proximity of the 
attractive glass, and at high temperature, because of the 
nearby repulsive glass. A power law divergence describes 
the attractive glass increase of rj with exponent 7^ ~ 3.75, 
i.e. the same that is found also for the density relaxation 
time 7t-. Data refer to an average over 20 independent 
starting configurations and over time for a minimum of 
200tq,. a pronounced reentrant behaviour, covering two 
full decades toward both limits, is observed in 77, similar 
to that reported previously for the diffusion coefficient D 
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FIG. 6: Viscosity approaching the attractive glass transition 
along path 2 (full black circles), and in the reentrant region 
along path 3 (empty red circles), as a function of attraction 
strength. Lines represent power-law fittings (with values of 
the critical attraction strength fixed to the previously deter- 
mined values reported in Table|Tl|, with exponents equal to 
3.16 for path 2 and 3.75 for the attractive side of the reentrant 
path 3. 



in the same systemf 




FIG. 7: ^ < (Ay4(t))^ > l%Vt for different attraction strength 
r along the isochore 0c = 0.58 for path 3A. On decreasing F, 
the long time limit first decreases (full lines) and then in- 
creases again (dashed lines), resulting in a pronounced reen- 
trant behaviour of the viscosity. 



V. COMPARISON OF Caaif) WITH MODE 
COUPLING THEORY 

MCT predicts [42j that the stress correlation function is 
related to an integral over all wavevectors of the density 
correlation functions: 



607r2 



dqq 



dlnSjq) 
dq 



(12) 



We theoretically calculate Ca-a{t) along two paths analo- 
gous to paths IB and 2 studied in simulations, to compare 
the full time-behaviour of the stress correlation function. 
Hence, we study: 

(i) a one-component hard sphere system with increasing 
(f), using the Percus-Yevick (PY) structure factor as in- 
put; 

(ii) a one-component AO model with size ratio g = 0.1 
at fixed packing fraction ip — 0.40. Here S{q) is calcu- 
lated using PY closure for the two-component Asakura- 
Oosawa mixture. This model mixture is composed of HS 
colloidal particles and ideal-gas polymers with HS inter- 
actions between polymers and colloids The obtained 
colloid-colloid structure factor is used as input to a one- 
component MCT, a treatment based on the validity of an 
effective one-component description for small polymer- 
colloid size ratio IHIHS]- We did not use the fundamental 
measure density functional theory [iGj |47] which yields 
analytical expressions for Sij{k) as done previously [4S] 
because within this closure the system shows spinodal 
instability before MCT would actually give a glass. This 
is not the case with PY closure for which only a very 
tiny increase in the structure factor at small q is found 
approaching the MCT transition. 

We solved the full dynamical MCT equations, as well 
as their long time limit, to calculate the viscoelastic prop- 
erties close to the glass transition. We used a grid a 1500 
wave- vectors with mesh Aq = 0.314. 

The long-time limit of the integrand of Eq. [12) 



t^oo dq 

(13) 

is plotted as a function of qa, in Figure|8]for both studied 
systems, being the critical non-ergodicity parameter at 
the MCT transition. The same figure reports also and 
the input static structure factor, also at the transition, 
S^{q). 

For the repulsive glass we find that the dominant con- 
tribution to the integral is provided by the wave-vector 
region around the nearest-neighbour peak, i.e. q*a ~ 6.5. 
For the attractive glass, on the other hand, the domi- 
nant contribution is found at much larger g- values, i.e. 
q* a fa 24: (in the region of the fourth peak of S{q)) pro- 
viding another confirmation of the importance of small 
length-scales in the localization properties of such a glass 
[¥!] . Moreover, in this case, the integrand is not just 
peaked around a specific value, but it is rather spread 



dlnSiq) 
dq 



rc 

J a 
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within a very large g-interval. The ampHtude of the in- 
tegrand is also much larger in the case of the attractive 
glass as compared to the repulsive glass. 



a small discrepancy at very long times. We attribute 
this difference to the fact that, in the case of attractive 
glasses, a large window of wavevectors contributes to the 
decay of the stress autocorrelation function (see Fig. [s]). 





FIG. 8: Mode coupling contributi ons to the viscosity 
/(g)/(607r^), with /(g) defined in Eq. 13 The wavevector 
at which I{q) is maximum, q*a, is « 6.5 for the repulsive 
glass and « 24 for the attractive glass. To compare, we re- 
port in the same figure also the g— dependence of the critical 
non-ergodicity parameter and of the static structure factor 
S^{q). 

We can then compare in the upper panel of Fig. |9]the 
theoretical stress correlation function with the squared 
theoretical density correlator 0^, (t) at the maximum of 
I{q)- We show two state points, one close to the re- 
pulsive glass and the other state close to the attractive 
one. Apart from an amplitude scaling factor, the dom- 
inant contribution is already sufficient to describe the 
long-time behaviour of Ca-a(t) for both attractive and re- 
pulsive glasses. However, for the attractive glass case, the 
decay of the squared density correlation shows a slightly 
smaller stretching as compared to Caa-{t), which causes 




FIG. 9: Stress correlation function Ccra{t) (full lines) for re- 
pulsive and attractive glasses calculated within MCT (top) 
and from simulations (bottom). Dashed lines are the squared 
density correlation functions (l>gt{t), arbitrarily scaled in am- 
plitude to overlap the long time behavior. For the MCT data, 
the wavevector q* is the one reported in Fig. [sj while in the 
simulation panel it is the one which provides the best long- 
time overlap between cj>q*{i) a-nd Crrait). 

In the lower panel of Fig. |9] the time dependence of 
both Caa{t) and (j)g*{t), as calculated from the simula- 
tion data, are also plotted. Here q* is the wavevector 
at which the agreement between the time dependence of 
Caa{t) and 0^. {t) is optimal. The q* values found in this 
way, respectively q* a « 7.5 and q*a ~ 26, agree very well 
with those predicted by the theory [30] ■ Moreover, the 
behaviour of Cacr{t) is well-described (within the numer- 
ical error) by a single squared density correlator for both 
glasses. The small discrepancy which was observed in 
the MCT data for the attractive glass is probably buried 
within the numerical noise. 

Finally we want to compare the elastic moduli for both 
glasses in the theoretical and numerical calculations. In 
order to calculate elastic and viscous moduli, the stress 
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C^a{0) A TO 


Tl P 


0P 


aa(0) 


A 




n p 


0.58 


181 0.18 0.024 


13.30 0.509 


0.42 


1650 


0.077 0.011 


81.48 0.325 


0.57 


156 0.16 0.026 


3.56 0.665 


0.41 


1506 


0.072 


0.011 


8.09 0.389 


0.55 


134 0.15 0.024 


1.18 0.759 


0.40 


1470 


0.061 


0.011 


3.49 0.585 


0.53 


83 0.23 0.025 


0.20 0.421 


0.39 


1404 


0.071 


0.012 


1.90 0.949 


0.50 


34 0.39 0.024 


0.03 0.353 


0.30 


724 


-0.085 


0.013 


0.07 1.757 



TABLE III: Parameters of the fitting of Caa{t) for states 
close to glass transition for soft-spheres (path IB). 



TABLE IV: Parameters of the fitting of Ccrcr(t) for states 
close to attractive glass transition (Path 2). 



correlation functions calculated from simulations have to 
be Fourier transformed: G{u!) = iojC{uj), where C{oj) is 
the Fourier transform of Cca{t). However, due to the 
noise in the correlation function, direct transformation 
produces very low quality results. Thus, we have fit- 
ted Cca{t) with empirical functional forms close to both 
glasses before performing the Fourier transform. We have 
chosen 



Path IB 
Path 2 



MCT 



181 32 
1650 127 



400 
6000 



3 

100 



TABLE V: Approximate values of initial value of the stress 
correlation value Ccrcr{0) and height of the plateau, /o- for 
paths IB and 2. The first two columns refer to simulation 
data and the last two to theoretical MCT predictions. 



cut) = a.(0){/(t/ro)+ 

A{l~f{t/To))eM-{t/rif}} 



(14) 



where f{x) is an even function that describes the short 
time relaxation of Ca-a{t)- f{x) = 1/(1 + x'^) for the 
repulsive glass (Fig. [2| and f{x) — exp{—x^} for the 
attractive glass (Fig. p| . tq represents a microscopic 
time scale, which should be state-independent, whereas 
Tl gives the time scale for the stress final relaxation. The 
parameter A gives the amplitude of the stored stress (so 
that ACcr„{0) is the height of the plateau in Ccra{t)) and (3 
is the stretching exponent, which according to the MCT 
prediction should be roughly equal to the stretching ex- 
ponent of the density-density correlation function at q* . 

In Table |III| we present the parameters of the fittings 
for Ccrcrit) for states along path IB, drawn in Fig. |2]as 
thin lines. As expected, tq is state-independent and ri 
increases substantially when the glass transition is ap- 
proached. A and /3 are correctly estimated only when 
the second relaxation is noticeable, i.e. above (/)c = 0.55; 
in these cases the amplitude is almost constant and /3 
is compatible with the value obtained from the density 
correlation function at g*, /3 = 0.52 [55] . 

The parameters of the fittings for the attractive glass 
(path 2), shown in Fig. [5] are given in Table IV As 
before, tq is almost constant, whereas ri increases dra- 
matically upon increasing the attraction strength. 

From the values of the fits, we can directly com- 
pare other quantities between theory and simulations: 
namely, the t — value of the stress correlation func- 
tion Ccrcr(O) and the height of the long-time plateau /o- 
for both glasses. The results from MCT and simula- 
tions are reported in Table |V] for both studied paths. 
For both glasses, the simulations provide a lower value 
of Cctct(O) and a larger value of fa with respect to MCT. 
Although numbers are not important per se when com- 
paring to MCT, the ratio fa/Cacr{0) is wrong by one 



order of magnitude for both attractive and repulsive 
glasses. This result seems to suggest that the factor- 
ization approximation [42] adopted to derive Eq 12 may 



be too severe, although the structural relaxation is ap- 
parently well described, as shown by the comparisons of 
Fig.d 

We finally directly compare the elastic and viscous 



moduli G'{lu) and G"{u}) in Fig. 10 for repulsive (top) 




10" 10 10 10 10 10" 10 10 10 10" 10 10" 10 

CO CO 




10"" 10 10 10 10" 10 10 10 10 10"" 10 10" 10 
CO CO 



FIG. 10: Shear moduli G' and G" from simulations (left) and 
MCT (right) for repulsive (top) and attractive glass (bottom). 
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and attractive glass (bottom). We observe qualitatively 
the same trends for both transitions in theory and simu- 
lations, despite a shift in the absolute numbers: 

(i) an increase of G'{uj) at large-cj (but smaller than 
the microscopic frequency) with the approach to the glass 
transition; 

(ii) the appearance of a minimum in G" which moves 
to lower and lower to with decreasing distance from the 
transition, in agreement with previous experimental and 
theroetical studies on both repulsive [121 HO] and at- 
tractive glasses |51l 152] . The minimum appears when 
e < 0.01 according to the theory (e — \Xg — X\/Xg, with 
X being either ^ or F) , and at slightly larger values of e 
according to the simulations; 

(iii) much larger moduli (up to one order of magnitude) 
for the attractive than for the repulsive glass. This obser- 
vation holds both for theory and simulations and agrees 
well with recent rheological measurements for thermo- 
reversible sticky spheres [19^ i53j. 

Overall, MCT correctly predicts the behavior of the 
viscoelastic properties on approaching both glass transi- 
tions. However, the results disagree again quantitatively, 
and more importantly in the ratio of the height of the 
plateau in G' (or minimum in G") with respect to G'^ 

(or G'^ax)- 



VI. BREAKDOWN OF STOKES-EINSTEIN 
RELATION 

Finally, we discuss the breakdown of the Stokes- 
Einstcin (SE) relation[5l [551 [Ml ISZl |5H1 [Ml ED] close 
to the glass transition for all different studied paths. 

We start by examining path I. Fig. [iT] shows the SE 
relation for the hard sphere binary system and the soft 
sphere poly disperse system. To allow for a unifying pic- 
ture, we plot the results as a function of the relative dis- 
tance to the estimated glass transition {(f>g — (j>). At low 
and moderate density, far from the transition the data 
are consistent with SE, although different values limits 
are obtained for model A or B; whereas the former takes 
the stick value, D-q/T = (STrtr)""'^, the latter goes to the 
slip limit: D-q/T = (27r cr)"^. The reason for this dif- 
ference is not clear [SH [HH US]. In both cases, as the 
system approached the glass transition, the SE relation 
breaks down significantly, both in the form D-q and Dt 
(see inset). 

Fig. [12| shows the SE relation for the attractive glass 
case (path II) and along the reentrance (path III). The 
former case is rather clean, and allows us to access a 
breakdown by two orders of magnitude with respect to 
the typical SE value, both in Drj/T and Dt (inset). For 
both paths, at large F (low T) a clear breakdown of both 
Dt and D-q/T is observed for the attractive glass. 

For path III (reentrance case) , one has to bear in mind 
that the path becomes parallel to the repulsive glass line 
at small F (see Fig. 1) and the increase is limited to the 
one observed in the HS case at the same packing. For 



Drj/T 




FIG. 11: Breakdown of the SE relation for Dq/T approaching 
the repulsive glass transition for paths lA (empty red circles) 
and IB (full black circles). For the hard sphere case, T = 1. 
Lines are guide to the eye. The two horizontal dashed lines 
mark the slip and stick values of the SE relation. Inset: Dt 
for the same paths. 



this path we have also performed BD simulations. The 
BD results, also shown in Fig. [12] coincide with the MD 
data at all state points investigated, confirming that the 
SE behavior close to both repulsive and attractive glass 
transitions does not depend on the microscopic dynamics. 



Data in Fig. 11 and Fig. 12 provide evidence that the 
breakdown of the SE is a phenomenon which can be ob- 
served in the vicinity of both the repulsive and the at- 
tractive glass transitions. Within the investigated state 
window, it appears that the magnitude of the breakdown 
is enhanced in the attractive glass case, speaking for 
the presence of more intense dynamical heterogeneities 
|(31l65l[55] when confinement is originated by short-range 
bonds rather than by the excluded volume caging. 



VII. CONCLUSIONS 

In this article we reported the behavior of the vis- 
cosity in two models for short-range attractive colloids 
along three different paths in the attraction-strength 
packing- fraction plane. Along the first path, the sys- 
tem approaches the repulsive hard-sphere glass transi- 
tion. Along the second path, it approaches the attractive 
glass. The third path is chosen in such a way that the 
system moves continuously from the repulsive to the at- 
tractive glass at constant packing fraction in the so-called 
re-entrant region [67j . In this case, we have also compared 
brownian and newtonian simulation results, confirming 
that the viscosity is independent on the microscopic dy- 
namics, in agreement with results based on the decay of 
density fluctuations in atomic liquids [68] . 

We find that the increase of the viscosity on approach- 
ing the glass transition is consistent with a power-law 
divergence. The divergence of -q can be described with 
the same exponent and critical packing fraction previ- 
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FIG. 12: Breakdown of the SE relation for Dr]/T approach- 
ing the attractive glass transition for paths 2 (circles) and 3 
(squares-MD and triangles-BD). Note the partial breakdown 
also at high T for the reentrant path due to the closeby repul- 
sive glass. The two horizontal lines mark the slip and stick 
SE values. Inset: Dt for the same paths. The star indicates 
the HS value for path 3. 

ously found for the collective relaxation time, but with 
an exponent different from the one that characterizes the 
divergence of the diffusion coefficient. This holds for both 
attractive and repulsive glass. 

As previously observed for diffusion and collective re- 
laxation, the viscosity shows a non monotonic behavior 
with the attraction strength in the reentrant region (path 
III), confirming once more the validity of the theoretical 
MCT predictions. 

To provide a connection between density relaxation 
and visco-elastic behavior we investigate the leading den- 
sity fluctuation contributions to the decay of the stress 
autocorrelation function within MCT. Interestingly, for 
the case of the repulsive glass, it is possible to identify a 
small range of wave-vectors (not far from the first peak of 
the structure factor) which are responsible for the visco- 
elastic behavior. In the case of the attractive glass, in- 
stead, the decay of the stress is associated to a much 
larger window of wave vectors, centered at much larger 
values. In this respect, the visco-elastic analysis con- 



firms that dynamic arrest is driven by the short-lengh 
scale introduced by the bonding. We also compare the 
simulation results for the frequency dependence of the 
elastic moduli with corresponding theoretical MCT pre- 
dictions, finding a substantial qualitative agreement. 

Finally, we have evaluated the Stokes-Einstein rela- 
tion. A clear breakdown of the relation is observed on 
approaching both glass lines, consistent with the differ- 
ent exponents characterizing the power-law dependence 
of diffusion and viscosity. The breakdown is particularly 
striking on approaching the attractive glass (a variation 
of the product D-q/T of up to two order of magnitude 
in the investigated range). Recent theoretical work on 
MCT seems to provide insights that could be useful to 
reconcile the decoupling of self-diffusion and viscosity (or 
relaxation time) within MCT[57j. It would be interesting 
in the future to deepen our knowledge of the connection 
between SE breakdown and the presence of dynamic het- 
erogeneities, which has been previously studied for the 
same model [51]. 

Note: While finalizing the manuscript, we become 
aware of a numerical study by Krekelberg et al. (cond- 
|mat/07050381) which also reports the non-monotonic be- 
havior of the viscosity along the reentrant path and the 
breakdown of the SE relation. In that work, Krekelberg 
et al. seek a connection between the structural and dy- 
namical properties of the system. We show here that 
MCT predicts correctly the properties of the system upon 
approaching the glass transitions, i.e. the connection be- 
tween structure and dynamics is the non-trivial one pro- 
vided by MCT. 
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